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Abstract 



o 
o 

We derive a family of singular iterated maps — closely related to Poincare maps — that describe chaotic 
interactions between colliding solitary waves. The chaotic behavior of such solitary wave collisions 
q depends on the transfer of energy to a secondary mode of oscillation, often an internal mode of the 

pulse. Unlike previous analyses, this map allows one to understand the interactions in the case when 
this mode is excited prior to the first collision. The map is derived using Mclnikov integrals and matched 
asymptotic expansions and generalizes a "multi-pulse" Melnikov integral and allows one to find not only 
multipulsc hctcroclinic orbits, but exotic periodic orbits. The family of maps derived exhibits singular 
behavior, including regions of infinite winding. This problem is shown to be a singular version of 
the conservative Ikeda map from laser physics and connections are made with problems from celestial 
mechanics and fluid mechanics. 

Solitary waves are solutions to time-dependent partial differential equations (PDE) 
which can be described as profiles of unchanging spatial shape which move at constant 
velocity. A fundamental question is what happens when two of them collide. In certain 
classes of PDE a phenomenon known as chaotic scattering is seen in collisions two waves 
will appear to bounce off each other several times and eventually move apart, with the 
number of "bounces" and the speed at which the waves eventually separate shown to be 
a very complicated function of the waves' initial speeds. This has been seen in numerical 
simulations dating back twenty-five years but never fully explained mathematically. We 

(f-) study a small system of ordinary differential equations, which may be derived as an ap- 

proximation to the PDE dynamics. This system has been shown to reproduce many of 
the features of the chaotic scattering and is much simpler to analyze. We show that the 

£ — . simplified system may be further reduced into an even simpler set of equations called an it- 

erated maps. These maps are then analyzed using the tools of dynamical systems. We find 
significant mathematical structure, including many bifurcations and "infinite horseshoes", 
in the iterated map that are responsible for chaotic scattering. 

1 Introduction 



Solitary waves — localized solutions to partial differential equations (PDEs) which translate at uniform 
velocity with a constant spatial profile — are a ubiquitous phenomenon in physical sciences. A fundamen- 
tal question relating to these objects is their behavior upon collisions, cither with other solitary waves, 
or with localized changes ("defects") to the medium through which they propagate. 

The dynamics of such collisions depends, of course, on the PDE in question, but there exist classes of 
equations in which qualitatively similar phenomena are expected. In strongly dissipative systems, two 
colliding solitary waves generally lose their distinct identities and merge into a single larger bump. At 
the opposite extreme are completely intcgrable systems whose solitary waves (solitons) survive collisions 
with their identities intact, due to a linear structure hidden deep in the underlying equations. This 
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Figure 1: (Color online) Kink-antikink collisions in Q. (a) 1-bounce solution at v = 0.27 > v c (b) 
capture at v = 0.21, (c) and (d) 2-bounce solutions at v = 0.1987 and v = 0.2268, (e) 3-bounce solution at 
v = 0.2228. The position of one kink, determined by fitting the numerical solution to a kink of undetermined 
position, is given by the dashed black lines. 



was first seen numerically by Zabusky and Kruskal |43j and confirmed by the discovery of the exact 
two-soliton solution [25], both for the Korteweg-de Vries equation. 

An interesting case where the spectrum of possible behaviors is much richer is that of systems which 
are both non-dissipative and non-integrable. Chaotic scattering between colliding solitary waves is a 
problem that has been rediscovered by numerous groups since first being hinted at numerical simulations 
in the 1970's pp. We defer a discussion of a history of the problem to the next section. 

We consider the particular example of the </> 4 equation which arises as a model problem in many 
areas of theoretical physics: 

4>u - <i>xx - 4> + 4> 3 = o. (i) 

This supports "kink" solitons of the form 

(f>{x,t) = ±tanh-5=, (2) 
v2 

where £ = x ^°~2 * and the kink velocity may take any value \v\ < 1, and the choice of a minus sign 
in the above formula describes the so-called antikink. Following [U [31 [5] 7 we numerically simulate the 
collision of a kink-antikink pair with equal and opposite velocity, and show the results in figure [T] 

For waves with initial speeds above some critical velocity v c the waves separate, albeit with reduced 
speed, figure [TJl. For most initial velicities below v c , the waves are captured and form into a single 
localized structure, figure [T|d. For initial velocities in certain "resonance windows," however, they may 
eventually separate, as in figure [TJ;, d, and e. Plots lc and Id show 'two bounce' solutions, as one can 
see from the dashed line that the kinks collide twice. One may discern that the transient bound state 
in figure [lfc oscillates twice (the brightest white spots in the middle) while the transient bound state 
in figure yjl oscillates three times. Figure [l^ shows a 3 bounce solution. In figure [2j we show how the 
escape velocity v ou t depends on the initial velocity vi n . The initial conditions corresponding to the 5 
experiments in figure 1 are marked in figure [2] Also plotted, using color, is the number of times the 
kink and antikink "bounce off" each other before escaping. In fact, we have found significant detail 
beyond what can be reasonably depicted in this figure. The structure is repeated at smaller widths with 
larger numbers of bounces before escape, and we have found numerical solutions with as many as nine 
bounces. Of course, as this PDE displays sensitive dependence on initial conditions, it also has sensitive 
dependence on numerical discretization and the fractal structure persists, but shifts slightly with small 
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Figure 2: (Color online) The output velocity as a function of the input velocity for kink-antikink collisions, 
color coded by number of times the kink and antikink "bounce" off each other before escaping. The labels 
(a)-(e) give the corresponding solutions in figure [Tj 



changes to the discretization. From figure [2] it is clear that the behavior of these colliding solitary waves 
is an example of chaotic scattering; see [33] for an introduction to a special issue dedicated to this topic. 

Notice in figure [I] that the escaping solution has an additional oscillation not present before the two 
solutions collide. In fact this oscillation is evidence of an internal mode which is excited during the 
collision, and which is responsible, as we will discuss below, for both the capture dynamics when v < v c 
and the subsequent escape in the case of resonance. Linearizing the 4 equation ( [Tj) abou t the kink ([2]), 
and looking for solutions of the form yi(£) cosier, where r = (t — v(x — xO))/vl — v 2 is time in the 
Lorenz-shifted coordinate frame, one finds solutions 

Xl(0 = ( ^ )1/2tanh ^l SeCh ^ (3) 

with frequency oj\ = \J\- A question of interest, then, is what happens if the kink's internal mode is 
excited before the collision. This is depicted in figure [3] The collision is initialized with 

(j>{x, 0) = cj> Q (x; x 0l 6*o, A)) ~ <Po(x; -xo, 0o, A ) - 1, 

where 

(f> (x; x , 9o, A) = <t>K (0 + Axi(£) cos (u x t - O ). 

The antikink and its internal mode are exact mirror images of the kink and its internal mode, and 
equation ([l]) preserves even symmetry. Each of the approximately 11500 pixels represents a numerical 
solution of 0. We did not consider the more complicated case where the kink, antikink, and their 
internal modes are chosen without this symmetry. The methods developed by Goodman and Haberman 
in [21 Q1)J [T5J [T7] are sufficient to describe much of the behavior depicted in figures [T] and [2] although 
the form of the map derived in the current paper for the more general class of initial conditions allows 
a much fuller understanding of the dynamics, especially of the structures seen in figure [3] 

The goal of this paper is to explain in detail the mechanisms underlying the chaotic scattering 
behavior described above and in particular to describe how to frame this system in such a way that the 
standard techniques of dynamical systems theory may be applied. We will proceed as follows. Section [2] 
contains further preliminaries including a historical overview of the phenomena, and an introduction of 
"collective coordinate" ordinary differential equation (ODE) systems used to model the PDE dynamics. 
In section [3] we derive, in detail, an iterated map that functions as a sort of Poincare map for the ODE 
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Figure 3: (Color online) For this figure a kink and antikink were initialized with incoming speed V[ n = 
0.19871, the center of the first two-bounce resonance in figure [2] and an internal mode pre-excited with 
parameters Aq and 0q. (a), the speed with which the kink and antikink separate; (b), the number of times 
the kink and antikink collide before escaping. 



model. We then interpret the map and show the results of direct numerical simulations in section [4] 
In section [5] we will analyze its dynamics. In section [6] we extend the analysis to two related systems, 
one from a different solitary- wave collision problem and the second from geophysical fluid mechanics 
and find interesting and more complex iterated maps. We conclude in section [7] with a discussion of 
how dissipation — caused in the PDE by the irreversible loss of energy to escaping radiation — may be 
introduced into the maps and used to explain capture. 



2 Further Preliminaries 

2.1 Background and technological motivation 

We have summarized the history of this problem in our previous work; see especially |17j . The two 
bounce resonance was discovered and explored in a series of papers by Campbell, Peyrard, and various 
collaborators in the 1980's [51 [H |2H |3B 29 following some hints in earlier numerical experiments by 
Ablowitz et al. pQ . Campbell et al. showed via careful examination of numerical experiments and via 
heuristic arguments that the capture and subsequent escape of the solitary waves was due to resonant 
energy transfer to a secondary mode, namely the kink's internal mode. Their calculations for the 
locations of the two-bounce windows, based on parameter-fitting to their numerical simulations, are 
remarkably accurate but do not give any insight into how to find v c or how the window locations might 
depend on parameters in the equations. The fractal structure was studied in somewhat more detail by 
Anninos et al [3]- 

The two-bounce resonance phenomenon was subsequently observed by Fei, Kivshar, and Vazquez in 
the collisions of kink- like solitary waves with localized defects [321 HI], and then by Tan and Yang in the 
collisions between solitary waves in a system of coupled nonlinear Schrodingcr equations describing light 
propagation in birefringent optical fibers |40l 142] . More recently it has been seen in quantum field theory 
in the interaction of topological solitons with defects in the metric of the background spacetime [2"51 155] , 
In section [6] we will extend our methods to the system studied by Fei et al. in [12] . 

There have been a few studies that looked at the behavior of solitary waves when the secondary 
mode is excited prior to the collision, for example in Fei et al. and Forinash et al [T5] . These have 
been mainly small-scale numerical simulations. One goal of this paper is to investigate such collisions 
more thoroughly. 

Many previous studies have used collective-coordinate models to study collision phenomena. In 
such models, the infinite-dimensional dynamics are reduced to a finite-dimensional system of ordinary 
differential equations (ODE). These ODE models are derived via the "variational method,' a non-rigorous 
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procedure based on the underlying Lagrangian structure of the system. Such methods are reviewed by 
Malomcd [28 . In recent work with Haberman, we have used these collective coordinate models to derive 
approximate formulas for the critical velocity and derived simple iterated maps that reproduce the fractal 
structure of figure [2] In section [3] we derive a more general form of this map which applies to a wider 
class of initial conditions and is more amenable to the methods of dynamical systems. 

In the earlier of the above-cited works, the PDE simulations typically used finite-difference (second or 
fourth order) discretization of the spatial derivatives and explicit Runge-Kutta time stepping. Figures II] 
and [2] were produced using pseudospectral (cosine-transform) spatial discretizations and time-stepping 
algorithms which treat the (stiff) linear portion of the equation exactly while applying an explicit fourth- 
order Runge-Kutta time step to the nonlinear terms [24 . This allows for the use of a coarser spatial 
mesh and much larger time steps, and, more importantly, eliminates the discretization-induced damping 
which may have made it much harder to resolve the narrow windows of chaotic scattering. 

One motivation that has often been cited for studying solitary wave collisions, especially collisions 
with defects, is the desire to build all-optical communications systems. Current technologies send signals 
at high speeds using pulses of light. Such signals are converted to electronic form for processing and 
then converted back to optical form for retransmission. The idea is that one may engineer the nonlinear 
effects in an optical medium to allow the data processing and thus avoid the time-consuming conversion 
back and forth to electronics. In one such scenario for an optical memory, a defect in an optical medium 
might be used to "trap" a pulse of light [2) [20]. At a later time, a second "probe" pulse might be sent 
in to detect whether a pulse has previously been captured. The computation depicted in figure [3] is 
meant to give an idea of the possible behavior of such a system. This may be seen as either a blessing 
or a curse: a curse because to ensure that a chaotic scatterer gives the desired output for a given input 
requires very accurate control on the input, and a blessing because, once obtaining such precision, there 
are a large number of possible outputs that could be used to encode different pieces of information. 



2.2 Collective Coordinate ODE Models 

Solutions to the <f> equation minimize the action with Lagrangian density 

m = \g - \€ + \^ - (4) 

A system of ODEs that models the behavior of the kink-antikink collision is derived using the "variational 
approximation." Instead of looking for minimizcrs of 

C(4>) dx dt, (5) 

we look for minimizers among functions of a predetermined spatial form dependent on a finite number of 
time-dependent parameters. For this particular problem, there are two such parameters: X(i) describing 
the separation between the kink and antikink, and A(t) giving the amplitude of the internal mode xi- 
This parameter-dependent profile is substituted into ^ , and the inner x integral is evaluated explicitly, 
giving an action dependent only on t. Computing the Frechet derivative of this action gives the Eulcr- 
Lagrange equations for the evolution of X(t) and A(t), 

X = 1 v ^ (-F(X)X 2 - U'(X) + F(X)'A) (6a) 



2 + 21 (X) 

A + ufA = F(X). (6b) 



We will study a simplification of these equations, so the exact forms of all the terms will not be 
important. The relevant details are summarized in figure Hi. The terms I(X) (position-dependence of 
mass) and F(X) (coupling) decay exponentially for large X\, while the potential U(X) has a single 
minimum and unbounded for X < and decays exponentially to U = 2 as X — > oo. Ignoring the 
coupling term F'(X)A in ( [6a|, the uncoupled X-dynamics conserves an energy E, whose level sets give 
the trajectories seen in figure |4b. This system was shown in simulations by Anninos et al. to reproduce 
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Figure 4: (Color online) (a) The functions denning system Q. (b) The phase-plane for equation (6a 
with the coupling term F'{X)A ignored. 



the chaotic scattering seen in the PDE system [3J and in [15], we showed how to determine the critical 
velocity and derived a restricted form of an iterated map that reproduces this behavior. 

Here, as in [T7], we simplify (|6| by eliminating I(X) and replacing U(X) and F(X) with simpler 
functions that preserve the main features of ([6]), namely the topology of the uncoupled phase plane. We 
have made this choice in order to clarify the exposition and will discuss later, where appropriate, how 
this effects the calculation. Our simplified model is 



mX + U'(X) + cF'(X)A = 
A + lu 2 A + eF(X) = 0, 



(7a) 
(7b) 



where 

U(X) = e~ 2X - e~ x and F(X) = e" x . 

U(X) is the potential for the Morse oscillator, a model for the electric potential of a diatomic model, used 
in quantum mechanics with exponentially weak long-distance interactions. We note that this system 
conserves an energy (Hamiltonian) of the form 



H = ™^X 2 + U(X) + i(i 2 + u?A 2 ) + eF{X)A. 



(8) 



That such an ODE can reproduce the qualitative and, to a surprising extent, the quantitative features 
of the kink-antikink collisions has been well-documented [3] [TTJ [TJ] 3DJ 02] ■ Figure [5ja) is the equivalent 
of figurc[2]for equation ([7]), with impressive qualitative agreement. A major difference between these two 
computations is that in the ODE model, capture only happens to a set of initial conditions of measure 
zero by conservation of phase-space volume, while in the PDE, a nonzero fraction of the initial conditions 
lead to capture. Each of the windows has a finite width, although most in this figure are unresolved. 
Figure |5^b) shows, as a function of both v and e, the number of interactions before the solution escapes 
to infinity. Figure [ga) of the fi gure may be thought of as a horizontal slice through Figure u5[b) . 

Figure [6] is the equivalent to figure [3] for a disk of initial conditions of equation Q . It is more 
informative — which we will show later on — to consider a disk of initial conditions of constant energy H 
given by ([8]). In this figure the point (x, y) — {A cos #0, A sin #0) corresponds to a numerical simulation 
to with 



X(0) - X n 



X(0) 



2H - uj 2 A 2 



A(0) = A cos (6>o - Widclay); 

A(0) = Ll>Aq Sin (6> - ^idolay), 



(9) 



where idciay is the time it takes a solution to (7al (neglecting A(t)) with (X(0),X(0)) as given to reach 
the manifold X — 0. The term tdeiay is included so that two adjacent initial conditions on a given radius 
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will arrive at the manifold X — with approximately the same phase, canceling a very rapid change of 
this phase as A — > \J2H /to, the outer edge of the diskj^] 

3 Derivation of the discrete map 

Perturbation methods are applicable when e -C 1. In addition to may be order-one or else may satisfy 
lu ^> 1. In model (TTu , both these conditions are met. Rigorous analysis in the case that e <C 1 and 
lu = 0(1) is standard and our results are essentially equivalent to those of Camassa et al. In the case 
that t < 1 and lu 3> 1, these results may be made rigorous by combining the above work with that 
of Delshams and Gutierrez [10]. In either case the formal calculation is the same. Our method is an 
extension of that used by Goodman and Haberman [l4 l \TE \ [TB I 117] . 

We let I = \{A 2 + lu 2 A 2 ) be the canonical action variable for the linear oscillator and consider two 
specific scalings. First, by setting e = in ([7]), the two equations decouple completely. 

Outer limit: mX + U'(X) = (10a) 

A + u 2 A = 0. (10b) 

We will apply this scaling when X ^> 1. The second scaling — the inner scaling — is to let A = eA'. In 
the limit e — > 0, this becomes 

Inner limit: mX + U'(X) = (11a) 

A' + lu 2 A' + F(X) = 0. (lib) 



The evolution of X(t) is identical in (10a I and ( 11a I . Its phase portrait contains a homoclinic 
orbit to a degenerate saddle point at X — oo which separates bounded, negative energy, solutions from 
unbounded, positive energy solutions, topologically equivalent to that shown in figure (4]}. A solution 
becomes trapped when e > and the coupling to the second mode A(t) causes a trajectory to cross from 
the region of unbounded trajectories to that of bounded, periodic, trajectories. 

Equation ^ conserves a time-dependent energy of the form Define E = ^X 2 + U(X),the 
energy in the mode X(t), so that the level set E = along the separatrix. The matched asymptotics 

1 This set of initial conditions does not lie precisely on a level set of the Hamiltonian, but one on which the H has variations 
of 0(e~ Xa,!lx ). This has very little effect on the figure. 
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Figure 6: (Color online) (a) The exit velocity and (b) the number of bounces before escape v ut as a 
function of the amplitude and phase of the secondary mode for initial condition ([9]). Run with parameters 
e = m = 1 and u = 2, and energy H = mVg/2, where vq = 0.287 is the velocity of the fourth two-bounce 
window in figure [5^ a). 



will depend on the assumption that the solution remains close to the separatrix throughout its evolution 
and will alternate between these two types of approximations. We define a sequence of times tj, which 
describes the instants at which the energy jumps between the two modes. For t — tj = O(l), the solution 

is well-approximated by the separatrix orbit, which has the explicit form Xs = log (l + ^r^~) > anc ^ 



2m 

that, backwards in time, the oscillatory mode of the approximate solution has the asymptotic behavior 

A(t) ^C(CjCoaw(t-tj) + Sjamu}(t-tj)), (12) 

where C is a constant to be defined momentarily The immediate objective is to find the asymptotic 
behavior of A(t) along the inner expansion, and to match this behavior, through the subsequent outer 
solution, to the behavior along the next inner solution, which is centered at a time tj+i, yet to be found. 
The linear oscillator can be solved, assuming the solution X can be approximated by the separatrix 
orbit X§, using the variation of parameters formula: 

A(t) =C(Cj cosui(t -tj) + Sj smuj(t -tj)) + 



e 



smu(t — tj) / F( Xq(t — tj)) cos lu(t — tj)d; 
U J-oo ( 13 ) 



e , . 
+ — coswlt — t, I 



I F(X s (t - tj)) sinw(r - tj)dr, 

J — OC 



so that as t — tj 



(14) 



A(t) - (cCj + ~ J F ( x s(r - tj)) sincj(r - t^dr^j cosuj(t - tj) 

+ (eSj - - J F(X s (t - tj)) cos cj(t - tj)dT^j sinuj(t - tj). 
As ^s(^) is an even function, the integral in the coefficient of cosu>(t — tj) vanishes identically. Letting 
C=-[ F(X s (t)) cos lut dr = - [ F(X s (r))e lWT dT, (15) 

^ J-oo U J-oo 
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then as t — t j — > oo , 

A(t) ~CC 3 cos Lj(t - tj) + C(Sj - 1) sin uj(t-tj). (16) 
Using the exact form of Xg and the residue theorem gives 

to J_ 00 2m + t z lu 
On the ensuing inner approximation, X(t) rs Xs(t — ij+i), and we write 

A{t) ~ C (Cj+i coscj(t — t, J+ i) + Sj+i smuj(t — tj+x)) , 



where tj+i remains to be found. We rewrite equation (16 1 in the limit as t — tj+± — > — oo using 

trigonometric identities and find = (_ C °L% + +! ™e£.) (sf-i) > where °o+i = ~ 

Letting Zj — Cj + iSj , this simplifies to 

Z j+1 = e- ie ^(Z j -i). (18) 

Before determining ij+i, we first examine the evolution of X(t) and its energy E(t). Along the 
jth outer solution (defined as occurring immediately preceding the jth inner solution), X(t) and A{t) 
are uncoupled, and Xj(t) lies along a level surface Ej = ^X^(t) + U(Xj(t)). Along the the inner 
solution that follows, we can compute dE/dt and thus AE, the change in energy between two successive 
approaches to the saddle point at infinity. We calculate 

^ = (mX + U'(X))X = -eAF'(X)X = -eA^-F(X(t)). 
dt dt 

Integrating this along the separatrix yields a Melnikov integral approximation to AE, since F(Xs(t)) — > 
as t — ► ±oo, 

/OO 7 />00 

A(t)-F(X s (t)) dt = e F(X s (t))A(t) dt. (19) 
-oo J — oo 



Using (13 1, this simplifies to 

AE= Up f(2S j -l). (20) 
Defining Ej = ^-£—£j, one finds (even before knowing 9j+i) that the map has a conserved Hamiltonian 



(derived from ([8])) 

H = S j + \Z j \ 2 . (21) 

An approximation to the interval (ty+i — tj) can be found using the matching condition for the inner 
and outer approximations to X(t). As t — tj — > oo, the separatrix satisfies 



x dY 



t-t j = ] /^J o -^L= = V2^(e x - 1)V» „ V2^ e ^ 2 + 0(e-^). (22) 

Along the near-saddle approach with energy Ej+\ and letting t = t* denote the time at which X assumes 
its maximum value X* , 



t* -t 



x ' dY 



2 J x y/E j+1 - U(Y) 
with the asymptotic expansion as t — t* — » — oo, 



' ' - ./-^i (^^ eW ) ~ " 2 ^ V ' 2 ' <23> 
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Summing (22 1 and (23 1 



t* -U 



2m tt 



—Ej + i 2 



A similar calculation yields an identical value for tj+\ — i*, so that combining them yields an asymptotic 
formula for the period 



9j+i - ujTj - Lo(t j+1 - tj) 
Thus, the full map may be written 



2m 



-E 



3+1 



2tt 
~C 



£j+i = Ej + 2 Imag Zj — 1 



Z 



3+1 



e -i9i+i( Zj -i). 



(24) 



(25) 



We may use conservation law (211 to eliminate £j from this map 



£ j+ i = £j + (2 Imag Z s -1) = H-Cj- Sj + 2Sj-l = H 
so the map may be written as a map from the complex plane to itself 



Z j+1 = eV^-» 2 -- H 



(Zj-i), 



where, from equations (17) and (24 1, a = = v2i 



2mw /e. Note that as either lj or e _1 will be 
assumed large, then a will also be large. Following [33], we change variables to Zj = Zj — | to produce 
our final form of the map: 



z J+1 ee HZ,) = e y/^-w-"( Zj - i/2) i/2. 



(26) 



Although this makes the map look slightly more complicated, this change of variables has the advantage 
that the map's inverse is of essentially the same form, 



Zj_ x =T-\Z j ) = eV'^ +l/2|2 - H (Z^/2)+z/2, 

and simplifies the formulas for fixed points that will follow. Defining the linear map p(Z) = Z* , the 
complex conjugate, then T and its inverse are related by: 

T~ x = p- l T P . (27) 

Remark 1 The approach taken here is in a sense orthogonal to the approach for two-dimensional 
Hamiltonian systems described by Guckenheimer and Holmes |21j and applied in our previous studies |181 
H9] . Writing the energy as H — Hq(X, X) + I + eHi(X, X , 9) where I and 9 are the canonical action- 
angle coordinates for the A-A coordinates, in that approach, one assumes that 4| > 0, and uses this 
fact to reduce the system by replacing the evolution variable t with 9. One then derives a Poincare map 
to the section 9 = #o- If, in the uncoupled case e = 0, the X-X system has a homoclinic orbit, one tries 
to show using a Melnikov integral that this homoclinic orbit persists in the Poincare map. One may 
consider T(Z) as a Poincare map, where the section is taken to be a subset of the hyperplane X = 0, 
switching the roles of X and A from the Guckenheimer and Holmes argument. 

Remark 2 The derivation of a similar map for the collective-coordinates model (J6j) proceeds in es- 
sentially the same manner as for our simplified model ([7]). Evaluation of the (Melnikov) integral (17) 
defining the scaling constant C is significantly harder and cannot be done in closed form, rather as an 
asymptotic expansion for large u>; see [T5]. Additionally, the approximation to the time interval contains 
an additional 0(1) at step p4|: 9 ^ 
term is absent from equation ( 24 1 . 



-£j+i + Ci for some constants C\ and C2- The C2-type 
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4 Interpretation and numerical iteration of the map 



4.1 Relation between map and ODE 



The map T defined by (26 1 greatly compresses the information content of the ODE system (J7j), and it 
is worth discussing the correspondence between the dynamics on the Z-plane and that of the full ODE. 



Map (26) depends on two constant parameters a and Ti. The total change of phase of a(t) between two 
consecutive interactions is directly proportional to a, and thus to lu. The constant Ti is a rescaling of the 
Hamiltonian H. Solutions on level sets Ti < may never escape, as the energy (|8| is positive-definite. 
Solutions on the level set Ti = may escape to oo, but only along the separatrix orbit, with escape 
velocity approaching zero. Almost every solution with Ti > escapes to oo at finite velocity as t — > ±oo. 



Map (26 1 is extremely similar to the Ikeda map, 

W j+1 = g(W j ) = be i 'U w iF+»+e. 

which arises in the modeling of lasers. Here < b < 1 is a dissipation factor, Wj £ C and e £ K. This 
is the composition of three operations: nonuniform rotation about the origin, contraction by a factor b, 
and translation by e. If b = 1, Q is the composition of two area- and orientation-preserving maps and 
thus preserves both area and orientation. 

Stolovitzky et al. studied a similar map of the form 

W j+1 = G{W j ) = be i '\ w ^ +e 

which, like the map T, is singular |39j . Map T is also formed as the composition of uniform transla- 
tion and nonuniform rotation — whether rotation takes place first, last, or in the middle, is equivalent 
dynamically. 

When Ti < 0, the map T is well-defined, continuous, and invertible on all of C. At Ti = 0, the 
rotation rate diverges at Z = i/2. When Ti > 0, the domain of map (26 1 is the complement of the disk 





i 




{• 













T is continuous outside 2? uti but the rotation rate diverges as the disk is approached. For points in 
2? out , the two-component form of the map, system (25 1, the next iterate £j+\ is well-defined and positive, 
but the angle 9 and thus likewise Cj+i and Sj+i 



are undefined. 
The map's range is the complement of the disk 



2? ln = 





i 




{" 













see figure [7] Note, from conservation law (21 1, when £ > 0, \Z + i/2\ 2 < Ti, meaning that points inside 
the disk T> ln correspond to portions of the trajectory outside the separatrix of figure [4Jd in differential 
equations ( |7a| . Similarly, when Z is outside the disk TD- m , the trajectory of (7a I are inside the separatrix. 
Note in figure [6|d, the set of initial conditions that escape on the first iterate (bounces = 1) is bounded 
by two arcs, the boundary of T> m and another that is nearly a circular arc on the interior of D- m . As 
e — * + , the approximation made in deriving T are become more accurate, and this second boundary 
becomes more circular. 

Of primary interest is the evolution of points in 2?i n . These points that correspond to solitons starting 
at a distance oo from the collision point. The map is iterated until the trajectory lands in the circle 
T> out and the soliton escapes. How many iterates this takes then tells how many bounces the solitary 
wave undergoes before escaping. The initial condition Z\ = —i/2 (at the center of 2?i n ) describes a 
trajectory along which A(t) — > as t — > — oo. If, for some n > 0, Z n = i/2 (at the center of X> ut) 
, then the solution escapes to X = oo after n in such manner that A(t) — > as t — > +oo. Thus, the 
energy level Ti contains an n-bounce resonance if f rn ~ 1 (—i/2) = i/2. Also of interest are the iterates 
initial conditions Z £ dV m . These correspond to trajectories along which (X, X) — » (oo, 0) as t — > — oo. 
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Figure 7: The disks V m and T> ont . When < TL < \, the disks do not intersect. When TL > \, the disks 
intersect. In (a) TL = 0.15, and in (b), TL = 0.5. 



Similarly if Z n E dT> out , the trajectory will approach (X, X) = (oo,0) as t — > +oo. The existence of a 
point Zi S <9£>i n such that J 7 ' 1-1 e dT> out , indicates that there exists an orbit homoclinic to a periodic 
orbit with energy TL. 

On the energy levels VTL > 1/2, the disks T>- ln and T> ont have an intersection of nonzero measure. 
Points in the intersection correspond to solitary waves of sufficient energy never to be trapped (they 
come from X — oo and return to X = oo with only one application of the map). The two points on 
&D m U dV out correspond to one-pulse homoclinic orbits to some periodic orbit. The energy level TL = 1 
is the energy level of the critical velocity for capture of an unexcited wave. At any higher energy the 
point Z = —i/2 which corresponds, recall, to the solitary wave arriving from X — oo with no energy in 
the internal mode is in T> out , and the soliton will escape without ever being captured. 



4.2 Numerical iteration 



We display numerical iterations of map J- . Figure [8] shows that map ( 26 1 can reproduce the figure [5^,, 
with parameters e = 0.25, m = 1, u) = 1. Each point on the a;-axis in figure |8[b) corresponds to the 
initial condition (£ ,Z ) = mvf n /oj 2 T^ 2 / , 0) in the map's two-component form ( |25| , while z; out oc V^n 
where the solution escapes to infinity on the n iteration. For these parameter, the critical velocity is 
computed correctly to about 9%, and the map calculation reproduces well the topological structure 
of the ODE simulations. Quantitative agreement between the map and the ODE simulations can be 

2 

improved by using the value of Ej + \ — Ej = — j 2 -, where v c is obtained from direct numerical simulation 
rather than from the Melnikov integral computation (19 1. 

Figure p\ shows the first four iterates of the disk T> ln under map (26 1, with the same parameter values 
as figure 18] and energy level corresponding to the initial velocity marked by ★ in that figure. This gives 
the parameter values a — 23.27 and TL — 0.324 in (26). Figure [l0|a) shows a portion of the curve 
T{d'Di ri ) which wraps around T> 



infinitely many times. This exhibits the stretching and folding typical 
of chaotic systems in addition to the singular behavior near near the outside edge of T> ln . Fgure 10 d, 



shows the number of iterates preceding escape depends sensitively on the initial condition in the disk 
2?i n , showing qualitative agreement with figure |6]x 
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Figure 9: (Color online) The initial conditions in the disk 2?; n , subfigure (a) and their first four iterates 
under map (|26|), subfigures (b)-(e). 
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5 Analysis of the iterated map 

5.1 Analogy to Sitnikov's reduced 3-body problem 

Moser considers the reduced 3-body problem, due originally to Sitnikov |32j . In this problem two pri- 
maries of identical mass M orbit about their joint center of mass in planar elliptical orbits of eccentricity 
e. A third body moves in a line normal to this plane and through the center of mass of the two primaries 
and evolves under attraction to the two primaries. The third body is assumed to exert negligible force 
on the first two. The third mass satisfies the differential equation 

(28) 



(z 2 + p 2 )3/ 2 ' 



where p = i — ecost + 0(e 2 ). 

He then considers the sequence of times t n at which z = and derives a map 



(Wl > v n+l ) = Se (t n > v n) > (29) 



where v n — |i(f n )|. In the limit e — > , the differential equation (28), and thus of the map Qq, as well, 
is completely intcgrable. In fact, the map Ga reduces to 



'Jn+l 



Vq; t n+1 = t n + T(v n ) =t Q + (n + l)T(v ), 



where the time delay T(v) is an increasing function of v such that T(0) = 0, T(v) — > oo as v — > 2 _ and 
T(v) is undefined for v > 2. Letting the variables t and v be polar coordinates for K 2 , Qq maps the circle 
of radius two to itself, with a twist angle that diverges as the circumference is approached. The map is 
undefined for v > 2, which correspond to orbits that travel monotonically from z — ±oo to z = 
and for which the notion of time between successive zeros is meaningless. When < e <C 1, the domain 
and range of Q e no longer coincide. The map Q t is defined on an ellipsoidal region Dq and maps onto 
a second ellipse D\ — GcDq. These regions are close analogs to the complements of T> out and T> ln . As 
e — > + , both regions approach the circle of radius two centered at the origin. 
Moser |32j proves the following theorem, after defining the integers 

tk+l — tk 
Sk = 9^ 
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Figure 11: The domain D$ and range D\ of map (29) 



as the number of complete revolutions made by the two primaries between two consecutive zeroes of 
z{t). 

Theorem 5.1 Given a sufficiently small eccentricity e > 0, there exists an integer m — m(e) such that 



any sequence s with s& > m corresponds to a solution to the differential equation (28 I 



Notice that the sequence Sk can be chosen completely arbitrarily. Further, one may define semi-infinite 
or finite sequences that begin or end with oo, corresponding to solutions that arrive from or escape to 
z = ±oo. This theorem is proven by constructing a horseshoe on which one defines a Bernoulli shift on 



a countably infinite number of symbols. In figure 12 we construct such a horseshoe for map (26 1 by 
considering the image of the topological rectangle 



jz : VR < ri < Z + % - < r 2 \ jz : VH < n < 





z-i 




< 




< r 2 j 




2 



Here we have set a = 15, n = 0.3, n = 1.2 
wraps about four times around the disk 2? Q ut 



% r 2 = 2Vn. The set F{1Z)) {^{K), respectively) 
(D in , respectively). T(JZ) intersects IZ in four stripes 



labeled Vi through Vg (H 2 through H 5 respectively). The subscript is given by [ot/yj \Z — i/2\ 2 — n\ 
for a point on the intersection of the stripe and the real Z axis and represents the number of complete 
oscillations made by A(t) between two consecutive bounces. 

The existence of this horseshoe implies the existence of a Cantor set of points which remain in the 
rectangle IZ for all iterations of T or T^ 1 on which the action of the map is equivalent to a Bernoulli 
shift on four symbols)^] As r\ / VH, the number of stripes contained in IZ grows without bound, so 
in this limit, the Bernoulli shift operator requires a countably infinite number of symbols. This is what 
allows one to specify the sequence of integers describing the number of complete oscillations between 
consecutive zeros in Theorem |5.1| Such "infinite horseshoes" are the subject of a recent preprint by 
Zambrano et al 

z ■ 



We may find an integrable limit of the map ( 26 1 by defining S 
which case 

Then as S 



l 



7 



H 



and y 



+ with 7 > fixed, map (26 1 approaches the integrable iteration 



In this limit T> out and D la degenerate to the unit circle. Letting X = y^ 1 , this map is identical to Gi 



o- 



V 2 



2 In fact to show the existence of chaos is sufficient to take the rectangle IZ smaller, so that it only contains the four stripes 
V3, H2 and H3. One may show the existence of a horseshoe without explicitly constructing it using Melnikov integral 



arguments alone, as in [9]. 
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Figure 12: The construction of the horseshoe, with the set J-(TZ) shown in light gray and the set J- 1 (1Z) 
in dark gray, (a) full image; (b) closeup on the rectangle 1Z. 



5.2 n-bounce resonant solutions and multipulse homoclinic orbits 

The analysis of section [3] can be used to determine some key features of the dynamics. Some of these 
results appear in our earlier work, but have a nice geometrical interpretation in the present context. 
First, the critical velocity for capture is found using (20 1. In the case that the internal mode is initially 

2 ,i2 

^-k—-, so that the critical velocity for capture is 



unexcited (Z = — i/2), one finds that AE - 



luC I ' \fm = eTr\/2e 



/2mu) 



We define an n-bounce resonant solution by the property Z — —i/2 and Z n _i = i/2, so that the map 
is undefined and the solitary wave escapes on the n iteration and no energy remains in the oscillating 
mode as t — > ±ooJ^] This may conveniently be written, using formula (27 1, as 



In the case n = 2m, 



and in the case n — 2m + 1 



jr-(m- 



= pF n 



(30) 



(31) 



(32) 



The two sides of equation (32 1 are complex conjugates of each other, so, equivalently 

T" 



6 



We illustrate this with explicit formulae for the 2- and 3-bounce initial velocities, and an implicit formula 
which yields the 4-bounce initial velocities. We first define the quantities (assuming the quantities under 



Note that the map as defined one component form (261 is undefined for Z n -\ = i/2, but that in the two component 
form — equations (24 1 and (25 I — this gives £ n = £q. 
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any square roots signs are positive) 



and 



<fi = a/yl — 7i, 



I 

2' 



lb = a/\ /4cos 2 H = a/ . /4cos 2 — . 

1 2 'V 2y/T^H 



- H. 



(33) 



Then the condition for a two-bounce resonant solution, from (31 1 is Z\ ■ 
solution is T{Z\) G M, and for four-bounce resonant solution is T{Z\) = Z\. 
Solving for the two-bounce resonance leads to the algebraic condition 



for the three-bounce resonance, 



and for the four-bounce resonance we find 



(2n- 1)tt; 



2n-l± 



> + ip= (2n- l)7r. 



for a three bounce resonant 

(34) 
(35) 
(36) 



Solving equations (34 1 and (35) for TL with fixed a allows us to find the two-bounce resonance 
velocities 



Auj 2 



V2,n 



(2n- l) 2 



and three-bounce resonance velocities 



«3,n± 



(37) 



(38) 



c (2n-l±|) 2 ' 

One cannot derive a closed-form expression for the four-bounce resonance velocities. 

Solutions to equation ( 30 1 on the energy level H = correspond to multi-pulse heteroclinic orbits of 
ODE system SolvingQ for a when H = gives a condition for the existence of a two-bounce 

homoclinic orbit to oo: 



a = V2u;e^ muJ = (2n - 1)tt 
and, using (35 1, the three-bounce homoclinic orbits satisfy 

a = V2ue^ rnuJ = \ 2n - 1 ± ij tt. 



(39) 



(40) 



Most interestingly, equation (36 1 becomes 

a |1 + 



1 

2 Icos ■ 



(2n- l)?r. 



The term 



diverges as a — > (2m — l)7r, which means that this equation has an infinite number 



2|cos f | 

of solutions as a — > (2m — l)7r from above or below. Thus, between any pair of two-bounce homoclinic 
orbits there exists a pair of three-bounce homoclinic orbits and a countably infinite sequence of four- 
bounce homoclinic orbits. This is summarized in figure [13] A similar pattern holds for two-, three-, and 
four-bounce resonances. 



4 To be more precise if such a solution exists, then, by the implicit function theorem, there exists a multipulse homoclinic 
orbit for a nearby value of a. 
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1 3 5 7 9 11 13 15 



Figure 13: The locations of two-, three-, and four-pulse homoclinic orbits, (squares, dots, and intersections 
between the solid curve and the dotted horizontal lines). 



5.3 Fixed points and periodic orbits 

Fixed points Zf of (26 1 must satisfy 



i.e. Zf € R. Using trigonometric identities and setting Zt — C;, we find the fixed points satisfy 



f 

1 (h 

■ - cot m 

2 V 2 



This has at most one root for each 



where 6i = a/ J Cf + \ — H, and we assume 2nl < 6i < 2ir(l + 1). 

value of I, and allows for a convenient indexing of the branches of fixed points Ci(a,H.). A more useful 
form for calculations is 

6 

2Ci sin ^ + cos ^ = 0. (41) 

The stability of a fixed point is determined by the eigenvalues of J{Zf), the Jacobian matrix^] 
Map (26 1 is orientation- and area-preserving, thus its Jacobian has unit determinant. Thus the stability 
is determined entirely by r — trace J{Zf). If |t| < 2, the fixed point is (neutrally) stable; if not, it 
is unstable. If r > 2, then Zf is a saddle point, and if r < —2, Zf is a saddle with reflection (both 
eigenvalues negative). At points where r = 2 one may find the canonical bifurcations, in this case 
saddle-nodes, and points where r = —2 correspond to period-doubling bifurcations. We may find most 
of these bifurcation values directly. 
The trace is given by 



a (Of 



,3/2 



(cf + l-n) 

At a saddle-node bifurcation, r = 2, which simplifies d42 



sin 6i + 2 cos 8i 



(42) 



to 



aR 



sin 



2 \2(R-H) 



_ cos __ sm _ 



(43) 



As map (|26[) is not an analytic function of 2. The Jacobian is taken with respect to the real and imaginary parts of Z. 
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The solution sin # = is inconsistent with (41 ). Setting the other factor zero is equivalent to 



n = cf + \ - ( a ( c? + - ) r, 



2/3 



At a flip (period-doubling) bifurcation, r — —2, which simplifies (42 1 to 



01 oiR . 9i 
cos — sin — 

2 I 2(R-H) 3/2 2 



(44) 



(45) 



This has solutions of two types. If cos y — 0, then = (21 + 1)tt and equation (41 1 further implies that 
Ci = 0. The locations of these bifurcations are then described by 



4 \(2l + l)n 



(46) 



Flip bifurcations can also happen where cos ^- / but the other term in (45 1 vanishes. This can be 
shown to imply 

1 f(Cf + l)a 



12 i lu \ 2 / 3 

H = Cf ' 1 ^ + lJ 



(47) 



4 V 4Ci 

At a flip bifurcation, a fixed point loses stability and a 2-cycle is created nearby. At bifurcations 



given by ( 47 1 , the flip bifurcation is supercritical, and a stable 2-cycle is created with zero imaginary 

(48) 



part and real part given by 




1 



4 (2/ + 1) 2 7T 2 



This 2-cycle is stable for 



1 



(2Z + 1)tt 



< H < - 
4 



4a 4 



(2Z + 1W (2Z + 1) 6 tt 6 ' 



(49) 



At this point, a stable 4-cycle appears. While we have not calculated further, we conjecture that a 
(Hamiltonian) period-doubling cascade occurs along each branch, leading to chaotic solutions as 7i is 
increased further. 

This is summarized in a bifurcation diagram showing the first few branches of Ci and C;.± for 
a = 20 and 7i varying, figure [14^, and in a bifurcation diagram over both parameters in figure [T4p . An 
infinite sequence of branches accumulates along the edge of the region Ti — C 2 — \ . The saddle node 
bifurcations approach the line C — as the branch index I — > oo, so that the range of parameters for 
which a fixed point is (neutrally) stable on a large-Z branch is very small. The 2-cycles which appear at 
the flip bifurcation described by ( |47| do not have a simple closed-form expression. Instead these 2-cycles 
T(Z + ) = Z- and T[ZJ) = Z + with Z_ = Z*, , the complex conjugate. Further Z± = C ±iS satisfies 



C 



— cot 



2Jc 2 + (s- \y -n 



-lir 



(50) 



c 



s- 



= — cot 



2Jc 2 + {s+\y-n 



— mn , 



where < l,m < it. A two-cycle resulting from the flip bifurcation (47 1 must, by continuity, satisfy 



I = m in ( 50 ) . Solutions to ( 50 1 may also arise in saddle-node bifurcations of this condition directly. In 



this case, they may satisfy I ^ m, as seen in figure 15 
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Figure 14: (Color online) (a) The first five branches on a bifurcation diagram for map (26 ) with a = 20. The 
bifurcations are marked: o, saddle- node (44); v> flip bifurcations of (46) A, flip bifurcations of (|47[) (these 
bifurcate into the complex iJ-plane, while all the branches shown lie on the real line). The black line 
is H. = C 2 + 1/4, the edge of the forbidden region 2? ut (in gray)- The marked points on the period-2 
curves indicate secondary period-doublings of (49). (b) A bifurcation diagram in both parameters, a and 
TC. The label 'Flip 1' refers to period-doubling bifurcations of type (46), and 'Flip 2' refers to those of 
type (47). Only the first five branches of saddle node and 'Flip 1' bifurcations are shown; subsequent 
branches approach the horizontal line TC = 1/4 for fixed a. 




Figure 15: (Color online) Bifurcation diagram for period-2 solutions of equation ( |50| ). (a), (b): the real 
and imaginary parts of three branches that arise from the period-doubling bifurcations marked with A in 
figure 14 with a = 20. (c), subfigd: the first 5 branches of period-two points to bifurcate due to saddle-node 
bifurcations of (50). A symmetric branch exists reflected over the real axis. The ordered pairs represent 



the numbers (l,m) in equation (50). 
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Comparison with ODE Bifurction Diagram 



We now compare the partial bifurcation diagram of figure 14 'a) with a partial bifurcation diagram of the 
underlying ODE ([7]). This will help us to interpret the map's dynamics and understand their domain of 
validity as an approximation to those of the ODE system. We first describe the relation between fixed 
points of the map p6| a nd periodic orbits of equation ([7]). Fixed points Zf = Ci of map (26 1 are shown 



preceding equation (411 to be real- valued. This implies, setting tj = in (|12|, that as t 
inner approximation, 

A(t) ~ C(Ci cos ojt + - sin tot) 



the 



Evaluating this approximation to A(t) at t = — T;/2, half of the period defined by equation (24 1, we find 



that 



= C 



cos 
sin 



2 



1)3, 



(51) 



Using the fixed-point equation (41 1 shows that — to this order in our approximation — d(— TJ/2) = 0. 
Since A(— 2}/2) = 0, this implies that the fixed points described by (41 1 correspond to solutions of Q 
even in t. Inverting the linear equation (51 1, 



(52) 



We numerically calculate periodic orbits of system ^ using a method of Viswanath that combines 
elements of Newton's method for root-finding and the Poincare-Lindstedt method for approximating 
closed orbits (41]. We use this method inside a pseudo-arclength continuation algorithm to compute 
how the period and shape of the periodic orbit changes as a function of energy level. Several branches 



of solutions are shown in figure 16 Subfigure 16 ^a) shows the value taken by a(t) on the periodic orbit 
where X reaches its maximum value (related to the fixed point on the section X = 0). There are seven 
branches, each plotted in a different color, five of which terminate in saddle-node bifurcations at small 
values of H. The two branches marked A and B, corresponding to the largest values of \a\, merge at 
H ps —0.258. Subfigure (b) shows the periods of those branches. The periods of branches A and B do not 
approach a common value — these branches do not merge in a simple saddle-node bifurcation. System [7] 
possesses a single fixed point (X*,a*) = (-log (2- 2u p-<? )- As iJ \ H(X*,a*) = -8/31, the 
periodic orbits along branches A and B shrink to the fixed point (X*,a*). The frequencies of these 
periodic orbits approach the two different frequencies associated with the linearization about (X*,a*). 
In this sense these two branches "end" at H = —8/31 this would not be possible in a bifurcation 
diagram describing simple fixed points. 

The a(t)-component of a periodic orbit of ODE Q is approximately given by 

a(t) « a c (t) + ao cos mt , 

where a c (t) is the complementary solution due to forcing from X(t) and has the same period as X(t). 
Along branch A, the cosine part of a(t) oscillates once for each oscillation of X(t), along branch B, twice. 
From subfigure (a), ao < along branch A, while ao > along branch B. Along the next branch CD, 
ao < 0. At point C, the cosine part of a(t) completes about two full oscillations, and as one moves along 
this branch toward point D an additional oscillation is added so that there are three full oscillations. At 
point E, the sign of ao flips from point D, and the number of oscillations grows from three to four as one 
move from point D to point E, this pattern continuing from branch to branch as the period is increased. 

Subfigures (c) and (d) show the quantity C; = a|x=x max coswXJ/2, defined in equation (52 1, which 
allows us to compare these calculations directly with figure 14 1. In subfigure (c), the branches are 



labeled using the same coloring scheme as in parts (a) and (b) . Part (d) is marked as in figure 14 1, with 
stability type calculated using the Floquet discriminant. This shows impressive agreement with that 
figure in some respects, but disagreement in others. In particular, the nearly "parabolic" branches are 
nested in the same manner in both figures, and each shows the same stability behavior — a small range 
of H for which the fixed point is neutrally stable. The "non-parabolic" branch corresponding to branch 
B is also quite similar, likewise stable (elliptic) far into the region of positive H, but the final branch in 
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Figure 16: (Color online) The seven branches of periodic orbits that cross the energy level H = 0. Colors 
are used in subfigures (a), (b), and (c) to distinguish branches, (a) a(t) evaluated at X = X max . (b) The 



period of each branch divided by 2tt. (c) Approximations to Cj, given by equation (52) to allow direct 



comparison to fixed points of map (26). (d) As in subfigure (c), but with the stability type indicated as 



in figure 14 1. 



the ODE figure has no analog for the discrete map. Further, no branch extends to H = — oo as is the 
case for the discrete map approximation. The derivation of map (26 1 depends on the solution staying 
close to the homoclinic orbit, a condition that is violated for large negative H as well as along periodic 
orbits near the fixed point — branches A and B. 



6 Extensions 

In this section, we consider extensions of the derivation of map ( 26 ) in section [3] The chaotic scattering 
phenomena described in sections [l] and [2] have been seen in a wide variety of systems and it is worth 



exploring the extent to which the explanation provided by map ( 26 1 and its analysis are sufficient to 
describe the chaotic scattering in these different settings, and to what extent the analysis needs to be 
modified. We here consider two related systems. 

The first is a model for kink-defect interactions in the sine-Gordon equation [T5J [TU Q~7] . It is nearly 
identical in form to Q above, with modifications to the potentials U and F, and explicit values for the 
constants: 

AX + U'(X) + F'(X)a = 0; (53a) 
a + X 2 a + eF(X) = (53b) 
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Figure 17: The unperturbed phase space, with the separatrix orbit in bold, for (a) the sine-Gordon 
model (53), and (b) the Lorenz model (54). In (b), the left and right edges are identified, forming a 
cylinder. Compare with figure |4^b) , which has only one separatrix. 



with potentials defined by 

U(X) 

and 



-2sech 2 (X) and F(X) = -2tanh(X) sech(X). 



X 2 = 



The last model we will consider was not derived in the context of solitary collisions, but is of essentially 
of the same form as the previous two. It was derived by Lorenz as a low-dimensional truncation of a 
set of "primitive equations" for atmospheric dynamics, and put in its present form by Camassa [27l 14], 
The equations take a very similar form: 



w 



-R sin ip, 



-z, z 



eR sin ip 



(54) 



where R is an order-one constant. This the canonical problem of a pendulum coupled to a linear spring. 
The difference between the three systems is topological, and the topological differences between the 



three systems lead to differing forms for the three maps. In both models (53) and (54), there exist two 
distinct separatrices as opposed to the single separatrix in model ([7]). The corresponding separatrix 
orbits Xg(t) are odd, and so are the coupling functions F(X). This asymmetry produces two slightly 
different maps depending on whether the solution is moving to the right or the left. In addition the the 
X — X phase space of the Lorenz model (54 1 is cylindrical. 



6.1 The sine-Gordon problem 



The method of deriving the map for problem (53 1 is identical to what has been described above. The 



one difference is that the phase space for the unperturbed X-dynamics features two separate homoclinic 
orbits. X — ±JCs(t) = ±sinh _1 (i), Because both Xg(t) and the coupling function F(X) are odd 
functions of their arguments, the maps obtained from the computation along the two orbits are different. 
If we assume that the kink is initialized traveling to the right from X = — oo, then if its initial velocity 
is below v c , it will follow the two heteroclinic orbits in alternation until eventually enough energy is 
returned to the propagating mode that it can escape back to X = ±oo. 



By properly scaling the variable a(t), one may show that system (53 1 has inner and outer scalings 



of essentially the same form as equations (10 1 and (11), allowing us to derive an iterated map in an 



analogous form to that in section [3j For solutions traveling along the upper separatrix, such that 
X(tj) = 0, and as t — > — oo, 



a(t) ~ C (Cj cosuj(t — tj) + Sj sinuj(t — tj)) 
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and C is again, for the moment, undetermined. Solving for a(t) by variation of parameters, as in 



equation (13 1 and allowing t — tj — >■ oo, one finds 



a(t) ~ [CCj + - / F(X s (t -t 3 )) sin X(t -t^dr ) cos X(t-tj) 



OO 
'OO 



CS'j - y / ^(^s(t - cos A(t - 4j)dT ) sin X(t - tj). 



— OO 



(55) 



A 

Because F(Xs(t) is odd, the integral in the coefficient of sin X(t — tj) vanishes. We define 



f C°° el" 00 
C = - / F(X s (t)) sin At dr = Imag - / F(X s (r))e iAr dr 

A J-00 A J-oo 



A ./-oo 1 + t 2 



using the exact form of Xg and the residue theorem. Then 

a(t) ~ C(C 3 + 1) cos A(i - tj) + CSj sin X(t -tj) as i - tj -> +oo. 

As in section [3] we recast this as the asymptotic behavior as t — tj+\ — > —oo for the next interaction. 
Using the same type of expansion as in that section, we find tj+i — tj ~ tt / yJ—2Ej + \. A calculation 
like (20 1 shows that AE = —2n 2 ee~ 2X (l — 2Cj). Letting, Ej = 2-K 2 ee~ 2X £j, the three dimensional map 



can be written 



£j+i — £j - (1 - 2C,-) 

C j+ i — cos8 j+ i(Cj - 1) + sin#,- + iSj (56) 
5j+i = — sin0j + i(Cj — 1) +cos^j + iS'j 



where 6* J+ i = A(tj+i — tj) = Xtt/ ^J—2Ej + i. Using the conserved energy H = £ j + Cj + Sj and again 
letting Zj = Cj + iSj yields the reduced map: 



= e J^^{Zj-l) (57) 



where a = Ae A /2y / e. Repeating this computation with the solution tracing along the leftgoing hete- 
roclinic orbit, the only change to map (57 1 is that wherever (Zj — 1) appears in the expression, it is 



replaced by (Zj + 1). 

If we assume that the solution first traces along a rightgoing heteroclinic orbit, then the general map 

is: 

Z j+1 = e Vl*i+c-° > j l 2 -* ( Zj + (-1)''). (58) 

This map has no fixed points, but one may look for solutions that satisfy Zj +2 = Zj and discover that 
this system does have fixed point, and one finds a bifurcation diagram very similar to figure [14^,. 

6.2 The Lorenz model 



Model (54 1 differs from the first two in two important ways. First, the evolution takes place on a cylinder. 
In the previous two models, the map became undefined if the solution to the ODE escaped to X — ±oo. 
In the present model, the analogous behavior is simply an orbit that wraps around the cylinder, in which 
case the map can be applied again. The second major difference is that the has genuine fixed points at 
ip = ±7r, of mixed elliptic-hyperbolic type in the full four-dimensional phase space, as opposed to the 
degenerate fixed points at infinity displayed by the first two models. This has a direct consequence on 
the general form of the asymptotic approximation of (tj+i — tj). 



Solutions to (54) conserve a Hamiltonian H = ^ R (cos ip + 1) + %- + The energy in the 

propagating mode is E = ^ — R 2 (cos ip + 1), where the constant has been chosen to make E = along 
the heteroclinic orbit. 
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The separatrix orbit in (54 1 is given by ips — ±sin ^anhitt. Assuming that the solution follows 
the upper heteroclinic orbit, centered at time tj, with 

x ~ C (Cj cos (t - tj) + Sj sin (t - tj)) 

then 

x =C(Cj cos (t - tj) + Sj sin (t - tj)+ 

eR 2 I — sin (t — tj) I sin ips( T — tj) cos (r — tj) dr + cos (t — tj) / sin ips(r — tj) sin (r — t/) dr I . 

V J —oc J— oc / 

As t — > oo, the integral in the coefficient of sin (i — tj-) vanishes. Once more using the exact form of 
and the residue theorem, we define 

C = eR 2 / sin (ip$ (r - tj)) sin (r - tj ) dr = eR 2 / 2 sech Rt tanh Rt sin t dt = 2-ire sech — 

where integration by parts and the residue theorem are used to evaluate the integral. Then as t — t j — > 

+00, 

x - C((Cj + 1) cos (t -tj) + Sj sm{t - tj)) . 
Now, assuming that the near-saddle approach occurs on the level set of Ej + i, we find 

2 , 4V2R 
tj+l ~tj = — log 



R *VWJ^\' 

This is defined regardless of the sign of -Ej+i, in contrast to the first two examples. This is a well known 
result describing the passage time near a saddle point and appears, for example, in Lichtenberg and 
Lieberman's textbook j^Hj. The rate at which At diverges as Ej — » is logarithmic, as opposed to the 
previous examples in which At cx E^ 1 / 2 . An equation like (19 1, yields 

E j+1 = Ej - 27r 2 esech 2 -^-(1 + 2Cj). 

2R 

Making the scaling Ej = 27r 2 esech 2 ^£j, defining Zj = Cj + iSj, and using the conserved quantity 
H = Ej + Cj + Sj , the map can be written as 



■ lop 



4V2R 



Z J+1 = e n \ c V\K-^+V 2 \J (Zj + 1) 

which is defined as long as \Zj + 1| 2 7^ H. 

If the solution traverses the lower separatrix orbit ip ~ —ips(t ~ tj), the above formula is modified 
by replacing (Zj + 1) with (Zj — 1). The maps are not defined in strict alternation, as for (58 1 — which 
formula applies is determined by whether the pendulum is swinging to the right or to the left at a given 
step. We can keep track of this by introducing a discrete variable Pj that can take only the values 
Pj = ±1. The state of the system can be described by the variables Zj and Pj where 



Z j+1 = e R \ C VI"-'W 2 I/ (Zj +Pj); (59a) 
P j+1 - sign (H - \Zj + Pj\ 2 )Pj. (59b) 

The rule for Pj says that if the energy £j+i = H — \Zj + Pj\ 2 is positive, then on the next pass, the 
pendulum will continue swinging in the same direction, but if it is negative, it will turn around and 
swing in the opposite direction on the next pass. 

This can be thought of as an iterated-map analog of a hybrid dynamical system. A formal definition 
of hybrid dynamical systems is given by Guckenheimer and Johnson |22j . Briefly, a hybrid dynamical 
system consists of an indexed collection of one or more ordinary differential equations. An initial 
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condition specifies both which equation is to be solved as well as what the initial values to be used for 
its solution. The solution is integrated until some time t\ when the trajectory crosses some specified 
manifold (if it ever does). At this point a discrete map is applied that tells the system which ODE 
to solve next and what initial conditions to use. This process then continues, producing a sequence of 
transition times tj and trajectories defined on intervals (tj-i,tj). Map (59 1 is the discrete time analog. 
At each step, one must determine the value of the parameter Pj in order to determine whether the 
pendulum is swinging clockwise or counterclockwise, and thus which form of the map to apply. The 
two maps are continuous (except on a singularity surface) maps, and a given map is iterated until it 
produces a negative value of £j+i, at which point, the evolution continues by iterating the other map. 
We do not know of any other examples of this type of map that have been studied. 

In the work of Camassa et al., the focus was on showing the existence of multipulsc homoclinic 
orbits |H|5]. Such homoclinic orbits take place on the energy level H — so that Pj+\ = —Pj and the 
novel aspects of the map are missing. The methods described in section 5.2 make finding multipulse ho- 
moclinics much easier. When TL > 0, the map becomes significantly more interesting, and equation (59) 
may be used to locate periodic orbits which may move clockwise and counterclockwise an arbitrary 
number of times in any sequence. 



7 Discussion and future directions 



In this paper we have studied an iterated-map model of solitary wave collisions. Our previous studies 
used a similar analytic method, but produced an iterated map that was significantly harder to analyze 
and which applied to only a limited set of initial conditions |14| I15[ 116} 117] . This paper analyzed in 
depth the map generated as an approximation to one such system. Chaotic scattering of a very similar 
type has been seen in many solitary wave collisions, and the methods used to analyze this map will not 
differ greatly, nor we believe, will the types of structures found, for example, in the bifurcation diagrams 
for such maps as we derived in section [6] 

We briefly discuss a few avenues not explored in the present study but of interest for future papers. 



Modeling of dissipation 

A key difference between the PDE and ODE simulations is the possibility of true capture for all time, 
as seen by comparing figures [2] and [5] By conservation of phase-space area, almost all solutions to 
equation Q that approach X = +oo as t — > — oo will eventually escape as t — > +oo. By contrast, the 
loss of kinetic energy to radiation acts as a damping to equation Q and makes capture possible for a 
significant set of initial conditions. 

One may modify equation ^ or ^ to include the effects of this dissipation. Such a dissipative 
model is derived perturbatively for the interaction of sine-Gordon kinks with localized defects in [15] , 



The effect was to modify the equivalent of equation ( 7b I to something similar to 

A + e 3 A 2 G(X)A + u?A + eF(X) = 0, 

where G(X) > and G(X) — > as X — * oo, i.e. the damping term is nonlinear in A and decays away 
from X = 0, so that damping only occurs when the kink and antikink are near each other. This should 



have the effect of modifying the equation for Ej in map (25) and destroying the conservation law (21 1. 
Thus, writing Hj = £j + \Zj\ 2 , we find Hj is strictly decreasing away from any fixed points. As Hj 
decreases, the escape regions V out shrink to a point and then disappear. If the disk disappears without 
the solution ever landing on it, then the solitary waves will be trapped forever. A detailed analysis of 
this dynamics is a topic for later research. Numerical simulations of the dissipative system in the above 
reference shows that this mechanism is enough to remove many but not all of the resonance windows. 
A systematic study of such modified maps is underway. 



Invariant manifold calculations 



In order to understand the fractal dynamics described in this paper, it is necessary to understand 
not merely the fixed points of map (26 1 but their stable and unstable manifolds, which control the 
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topological organization of the phase space. The subject of lobe dynamics, first developed by Rom- 
Kedar to understand the topology of such stable manifolds [37J HE] have the potential to be of great 
use. First of all, using measure-theoretic extensions of these methods due to Meiss, should be useful 
for calculating quantities such as the average number of times two solitons collide before escaping [30 . 
Recent work by Mitchell and Delos [31 allows one to gain a detailed understanding of the topology of 
the invariant manifolds and may be useful for gaining a more complete picture of the interleaving of the 
different scattering trajectories, far beyond the simple Cantor-set descriptions of section |5~T| 

This analysis is complicated in two ways by the infinite winding around the disks T> m and T> ont in 
map (26 1 with TL > 0. First, it implies the existence of an infinite number of fixed points, whose invariant 
manifolds may be important to the complete application of the above methods. Secondly, the invariant 
manifolds which cross these disks are split into an infinite number of disconnected pieces, so standard 
methods, which assume that these manifolds are continuous, are not directly applicable. It may be 
possible to avoid this difficulty by using different types of iterated-map reductions than equation ( 26 1 , 
or it may be necessary to modify the methods in order to make them directly applicable to this map. 



Application of methods directly to PDE simulations 

One of the applications often discussed for solitary wave-defect collisions is to build optical components 
such as switches or logic gates with them. Without getting into details, it would be necessary to 
understand figures such as figure [3] for those particular systems in order to tune initial conditions to be 
sufficiently far away from the boundaries between the scattering regions. To accomplish this, it should 
be necessary to implement many of the ideas of this paper directly to Poincare maps built directly from 
numerical simulations of partial differential equations. 
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